### Plot: Overview hybrid models

rm(list=ls(all=TRUE))


#packages
library(tidyverse)
library(sjmisc)
library(sjPlot)
library(sjlabelled)
library(ggplot2)
library(gghighlight)
library(RColorBrewer)
library(cowplot)
library(data.table)
library(dotwhisker)
library(ggpubr)
library(gridExtra)
library(lemon)
library(grid)

select <- dplyr::select
filter <- dplyr::filter


##### HYBRID MODELS ####

##### read data ####
setwd("SET_YOUR_PATH/")

#raw hhinc
bhps_r <- read_stata("./BHPS/bhps_coef_hhinc.dta")
gles_r <- read_stata("./GLES/gles_coef_hhinc.dta")
gss_r <- read_stata("./GSS/gss_coef_hhinc.dta")
liss_r <- read_stata("./LISS/liss_coef_hhinc.dta")
nlsy97_r <- read_stata("./NLSY97/nlsy97_coef_hhinc.dta")
polat_r <- read_stata("./POLAT/polat_coef_inc.dta")
psid_r <- read_stata("./PSID/psid_coef_hhinc.dta")
shp_r <- read_stata("./shp/shp_coef_hhinc.dta")
soep_r <- read_stata("./SOEP/soep_coef_hhinc.dta")

#hhinc_dec_pp_sqrt
bhps <- read_stata("./BHPS/bhps_coef_hhinc_dec_pp_sqrt.dta")
gles <- read_stata("./GLES/gles_coef_hhinc_dec_pp_sqrt.dta")
gss <- read_stata("./GSS/gss_coef_hhinc_dec_pp_sqrt.dta")
liss <- read_stata("./LISS/liss_coef_hhinc_dec_pp_sqrt.dta")
nlsy97 <- read_stata("./NLSY97/nlsy97_coef_hhinc_dec_pp_sqrt.dta")
polat <- read_stata("./POLAT/polat_coef_inc_dec.dta")
psid <- read_stata("./PSID/psid_coef_hhinc_dec_pp_sqrt.dta")
shp <- read_stata("./shp/shp_coef_hhinc_dec_pp_sqrt.dta")
soep <- read_stata("./SOEP/soep_coef_hhinc_dec_pp_sqrt.dta")

#loginc
bhps_ln <- read_stata("./BHPS/bhps_coef_loginc.dta")
gles_ln <- read_stata("./GLES/gles_coef_loginc.dta")
gss_ln <- read_stata("./GSS/gss_coef_loginc.dta")
liss_ln <- read_stata("./LISS/liss_coef_loginc.dta")
nlsy97_ln <- read_stata("./NLSY97/nlsy97_coef_loginc.dta")
polat_ln <- read_stata("./POLAT/polat_coef_loginc.dta")
psid_ln <- read_stata("./PSID/psid_coef_loginc.dta")
shp_ln <- read_stata("./shp/shp_coef_loginc.dta")
soep_ln <- read_stata("./SOEP/soep_coef_loginc.dta")

#hhinc_sat
bhps_subj <- read_stata("./BHPS/bhps_coef_hhinc_sat.dta")
gles_subj <- read_stata("./GLES/gles_coef_hhinc_sat.dta")
gss_subj <- read_stata("./GSS/gss_coef_sat_fin.dta")
liss_subj <- read_stata("./LISS/liss_coef_hhinc_sat.dta")
polat_subj <- read_stata("./POLAT/polat_coef_satinc.dta")
shp_subj <- read_stata("./shp/shp_coef_hhinc_sat.dta")
soep_subj <- read_stata("./SOEP/soep_coef_hhinc_sat.dta")


####p repare data ####

# delete empty rows
bhps_r <- bhps_r[-13,]
gles_r <- gles_r[-23,]
gss_r <- gss_r[-5,]
liss_r <- liss_r[-11,]
nlsy97_r <- nlsy97_r[-5,]
polat_r <- polat_r[-19,]
psid_r <- psid_r[-3,]
shp_r <- shp_r[-7,]
soep_r <- soep_r[-11,]

bhps <- bhps[-13,]
gles <- gles[-23,]
gss <- gss[-5,]
liss <- liss[-11,]
nlsy97 <- nlsy97[-5,]
polat <- polat[-19,]
psid <- psid[-3,]
shp <- shp[-7,]
soep <- soep[-11,]

bhps_ln <- bhps_ln[-13,]
gles_ln <- gles_ln[-23,]
gss_ln <- gss_ln[-5,]
liss_ln <- liss_ln[-11,]
nlsy97_ln <- nlsy97_ln[-5,]
polat_ln <- polat_ln[-19,]
psid_ln <- psid_ln[-3,]
shp_ln <- shp_ln[-7,]
soep_ln <- soep_ln[-11,]

bhps_subj <- bhps_subj[-13,]
gles_subj <- gles_subj[-23,]
gss_subj <- gss_subj[-5,]
liss_subj <- liss_subj[-11,]
polat_subj <- polat_subj[-19,]
shp_subj <- shp_subj[-7,]
soep_subj <- soep_subj[-11,]


#set new colnames 
colnames(bhps_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(gles_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(gss_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(liss_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(nlsy97_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(polat_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(psid_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(shp_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "id")
colnames(soep_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")

colnames(bhps) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(gles) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(gss) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(liss) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(nlsy97) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(polat) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(psid) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(shp) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "id")
colnames(soep) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")

colnames(bhps_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(gles_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(gss_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(liss_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(nlsy97_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(polat_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(psid_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(shp_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "id")
colnames(soep_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")

colnames(bhps_subj) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(gles_subj) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(gss_subj) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(liss_subj) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(polat_subj) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")
colnames(shp_subj) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "id")
colnames(soep_subj) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N",  "id")



#add identifiers for datasets
bhps_r$df <- 1
gles_r$df <- 2
gss_r$df <- 3
liss_r$df <- 4
nlsy97_r$df <- 5
polat_r$df <- 6
psid_r$df <- 7
shp_r$df <- 8
soep_r$df <- 9

bhps$df <- 1
gles$df <- 2
gss$df <- 3
liss$df <- 4
nlsy97$df <- 5
polat$df <- 6
psid$df <- 7
shp$df <- 8
soep$df <- 9

bhps_ln$df <- 1
gles_ln$df <- 2
gss_ln$df <- 3
liss_ln$df <- 4
nlsy97_ln$df <- 5
polat_ln$df <- 6
psid_ln$df <- 7
shp_ln$df <- 8
soep_ln$df <- 9

bhps_subj$df <- 1
gles_subj$df <- 2
gss_subj$df <- 3
liss_subj$df <- 4
polat_subj$df <- 6
shp_subj$df <- 8
soep_subj$df <- 9

#add model identifiers
bhps_r$model <- rep(1:2, times=6, each=1)
gles_r$model <- rep(1:2, times=11, each=1)
gss_r$model <- rep(1:2, times=2, each=1)
liss_r$model <- rep(1:2, times=5, each=1)
nlsy97_r$model <- rep(1:2, times=2, each=1)
polat_r$model <- rep(1:2, times=9, each=1)
psid_r$model <- rep(1:2, times=1, each=1)
shp_r$model <- rep(1:2, times=3, each=1)
soep_r$model <- rep(1:2, times=5, each=1)

bhps$model <- rep(1:2, times=6, each=1)
gles$model <- rep(1:2, times=11, each=1)
gss$model <- rep(1:2, times=2, each=1)
liss$model <- rep(1:2, times=5, each=1)
nlsy97$model <- rep(1:2, times=2, each=1)
polat$model <- rep(1:2, times=9, each=1)
psid$model <- rep(1:2, times=1, each=1)
shp$model <- rep(1:2, times=3, each=1)
soep$model <- rep(1:2, times=5, each=1)

bhps_ln$model <- rep(1:2, times=6, each=1)
gles_ln$model <- rep(1:2, times=11, each=1)
gss_ln$model <- rep(1:2, times=2, each=1)
liss_ln$model <- rep(1:2, times=5, each=1)
nlsy97_ln$model <- rep(1:2, times=2, each=1)
polat_ln$model <- rep(1:2, times=9, each=1)
psid_ln$model <- rep(1:2, times=1, each=1)
shp_ln$model <- rep(1:2, times=3, each=1)
soep_ln$model <- rep(1:2, times=5, each=1)

bhps_subj$model <- rep(1:2, times=6, each=1)
gles_subj$model <- rep(1:2, times=11, each=1)
gss_subj$model <- rep(1:2, times=2, each=1)
liss_subj$model <- rep(1:2, times=5, each=1)
polat_subj$model <- rep(1:2, times=9, each=1)
shp_subj$model <- rep(1:2, times=3, each=1)
soep_subj$model <- rep(1:2, times=5, each=1)


###recode DVs

#raw hhinc
bhps_r$dv <- rec(bhps_r$id, rec = "1=3 [Political Interest]; 2=1 [Vote]; 4=2 [Intention to vote]; 
               3=4 [Internal Efficacy]; 5=5 [Duty to vote];  6=6 [Vote: pers. benefits]")
bhps_r$term <- as_character(bhps_r$dv)

gles_r$dv <- rec(gles_r$id, rec = "1=1 [Vote]; 2=3 [Political Interest]; 3=2 [Intention to vote]; 
               4=4 [Has Party ID]; 5=5 [Party ID strength]; 6=10 [Internal Efficacy]; 7=11 [Imp.: Elections]; 
               8=12 [Imp.: Elec. Campaign]; 9=13 [Pol. media use]; 
               10=15 [Pol. Know.: Verbal]; 11=16 [Pol. Know: Visual]")
gles_r$term <- as_character(gles_r$dv)

gss_r$dv <- rec(gss_r$id, rec = "1=29 [GSS: Interest intl. affairs]; 2=30 [GSS: Interest mil. policy]")
gss_r$term <- as_character(gss_r$dv)

liss_r$dv <- rec(liss_r$id, rec = "1=3 [Political Interest]; 2=1 [Vote]; 3=2 [Intention to vote];  
               4=4 [Internal Efficacy]; 5=5 [Pol. part. index]")
liss_r$term <- as_character(liss_r$dv)

nlsy97_r$dv <- rec(nlsy97_r$id, rec = "2=2 [NLSY97: Political Interest]; 1=1 [NLSY97: Vote]")
nlsy97_r$term <- as_character(nlsy97_r$dv)

polat_r$dv <- rec(polat_r$id, rec = "3=2 [Political Interest]; 1=1 [Intention to Vote]; 2=7 [Pol. media use]; 4=3 [Internal Efficacy]; 
                5=6 [Duty to vote]; 6=4 [Has party ID]; 7=5 [Party ID strength]; 8=8 [Pol. part. index];  9=9 [Pol. part. onl.]")
polat_r$term <- as_character(polat_r$dv)

psid_r$dv <- rec(psid_r$id, rec = "1=1 [PSID: Vote]")
psid_r$term <- as_character(psid_r$dv)


shp_r$dv <- rec(shp_r$id, rec = "1=1 [Vote]; 2=2 [Political Interest]; 3=3 [No. polls part.]")
shp_r$term <- as_character(shp_r$dv)

soep_r$dv <- rec(soep_r$id, rec = "1=1 [Vote]; 2=3 [Political Interest]; 3=25 [Has party ID]; 
               4=2 [Intention to vote]; 5=26 [Party ID strength]")
soep_r$term <- as_character(soep_r$dv)


us_r <- bind_rows(gss_r, nlsy97_r, psid_r)
d_r <- bind_rows(bhps_r, gles_r, gss_r, liss_r, nlsy97_r, polat_r, psid_r, shp_r, soep_r)



#hhinc deciles (eq)
bhps$dv <- rec(bhps$id, rec = "1=3 [Political Interest]; 2=1 [Vote]; 4=2 [Intention to vote]; 
               3=4 [Internal Efficacy]; 5=5 [Duty to vote]; 6=6 [Vote: pers. benefits]")
bhps$term <- as_character(bhps$dv)

gles$dv <- rec(gles$id, rec = "1=1 [Vote]; 2=3 [Political Interest]; 3=2 [Intention to vote]; 
               4=4 [Has Party ID]; 5=5 [Party ID strength]; 6=10 [Internal Efficacy]; 7=11 [Imp.: Elections]; 
               8=12 [Imp.: Elec. Campaign]; 9=13 [Pol. media use]; 
               10=15 [Pol. Know.: Verbal]; 11=16 [Pol. Know: Visual]")
gles$term <- as_character(gles$dv)

gss$dv <- rec(gss$id, rec = "1=29 [GSS: Interest intl. affairs]; 2=30 [GSS: Interest mil. policy]")
gss$term <- as_character(gss$dv)

liss$dv <- rec(liss$id, rec = "1=3 [Political Interest]; 2=1 [Vote]; 3=2 [Intention to vote];  
               4=4 [Internal Efficacy]; 5=5 [Pol. part. index]")
liss$term <- as_character(liss$dv)

nlsy97$dv <- rec(nlsy97$id, rec = "2=2 [NLSY97: Political Interest]; 1=1 [NLSY97: Vote]")
nlsy97$term <- as_character(nlsy97$dv)

polat$dv <- rec(polat$id, rec = "3=2 [Political Interest]; 1=1 [Intention to Vote]; 2=7 [Pol. media use]; 4=3 [Internal Efficacy]; 
                5=6 [Duty to vote]; 6=4 [Has party ID]; 7=5 [Party ID strength]; 8=8 [Pol. part. index];  9=9 [Pol. part. onl.]")
polat$term <- as_character(polat$dv)

psid$dv <- rec(psid$id, rec = "1=1 [PSID: Vote]")
psid$term <- as_character(psid$dv)


shp$dv <- rec(shp$id, rec = "1=1 [Vote]; 2=2 [Political Interest]; 3=3 [No. polls part.]")
shp$term <- as_character(shp$dv)

soep$dv <- rec(soep$id, rec = "1=1 [Vote]; 2=3 [Political Interest]; 3=25 [Has party ID]; 
               4=2 [Intention to vote]; 5=26 [Party ID strength]")
soep$term <- as_character(soep$dv)


us <- bind_rows(gss, nlsy97, psid)
d <- bind_rows(bhps, gles, gss, liss, nlsy97, polat, psid, shp, soep)



#log income (eq)
bhps_ln$dv <- rec(bhps_ln$id, rec = "1=3 [Political Interest]; 2=1 [Vote]; 4=2 [Intention to vote]; 
               3=4 [Internal Efficacy]; 5=5 [Duty to vote]; 6=6 [Vote: pers. benefits]")
bhps_ln$term <- as_character(bhps_ln$dv)

gles_ln$dv <- rec(gles_ln$id, rec = "1=1 [Vote]; 2=3 [Political Interest]; 3=2 [Intention to vote]; 
               4=4 [Has Party ID]; 5=5 [Party ID strength]; 6=10 [Internal Efficacy]; 7=11 [Imp.: Elections]; 
               8=12 [Imp.: Elec. Campaign]; 9=13 [Pol. media use]; 
               10=15 [Pol. Know.: Verbal]; 11=16 [Pol. Know: Visual]")
gles_ln$term <- as_character(gles_ln$dv)

gss_ln$dv <- rec(gss_ln$id, rec = "1=29 [GSS: Interest intl. affairs]; 2=30 [GSS: Interest mil. policy]")
gss_ln$term <- as_character(gss_ln$dv)

liss_ln$dv <- rec(liss_ln$id, rec = "1=3 [Political Interest]; 2=1 [Vote]; 3=2 [Intention to vote];  
               4=4 [Internal Efficacy]; 5=5 [Pol. part. index]")
liss_ln$term <- as_character(liss_ln$dv)

nlsy97_ln$dv <- rec(nlsy97_ln$id, rec = "2=2 [NLSY97: Political Interest]; 1=1 [NLSY97: Vote]")
nlsy97_ln$term <- as_character(nlsy97_ln$dv)

polat_ln$dv <- rec(polat_ln$id, rec = "3=2 [Political Interest]; 1=1 [Intention to Vote]; 2=7 [Pol. media use]; 4=3 [Internal Efficacy]; 
                5=6 [Duty to vote]; 6=4 [Has party ID]; 7=5 [Party ID strength]; 8=8 [Pol. part. index];  9=9 [Pol. part. onl.]")
polat_ln$term <- as_character(polat_ln$dv)

psid_ln$dv <- rec(psid_ln$id, rec = "1=1 [PSID: Vote]")
psid_ln$term <- as_character(psid_ln$dv)

shp_ln$dv <- rec(shp_ln$id, rec = "1=1 [Vote]; 2=2 [Political Interest]; 3=3 [No. polls part.]")
shp_ln$term <- as_character(shp_ln$dv)

soep_ln$dv <- rec(soep_ln$id, rec = "1=1 [Vote]; 2=3 [Political Interest]; 3=25 [Has party ID]; 
               4=2 [Intention to vote]; 5=26 [Party ID strength]")
soep_ln$term <- as_character(soep_ln$dv)


us_ln <- bind_rows(gss_ln, nlsy97_ln, psid_ln)
d_ln  <- bind_rows(bhps_ln, gles_ln, gss_ln, liss_ln, nlsy97_ln, polat_ln, psid_ln, shp_ln, soep_ln)



#subjective
bhps_subj$dv <- rec(bhps_subj$id, rec = "1=3 [Political Interest]; 2=1 [Vote]; 4=2 [Intention to vote]; 
               3=4 [Internal Efficacy]; 5=5 [Duty to vote]; 6=6 [Vote: pers. benefits]")
bhps_subj$term <- as_character(bhps_subj$dv)

gles_subj$dv <- rec(gles_subj$id, rec = "1=1 [Vote]; 2=3 [Political Interest]; 3=2 [Intention to vote]; 
               4=4 [Has Party ID]; 5=5 [Party ID strength]; 6=10 [Internal Efficacy]; 7=11 [Imp.: Elections]; 
               8=12 [Imp.: Elec. Campaign]; 9=13 [Pol. media use]; 
               10=15 [Pol. Know.: Verbal]; 11=16 [Pol. Know: Visual]")
gles_subj$term <- as_character(gles_subj$dv)

gss_subj$dv <- rec(gss_subj$id, rec = "1=29 [GSS: Interest intl. affairs]; 2=30 [GSS: Interest mil. policy]")
gss_subj$term <- as_character(gss_subj$dv)

liss_subj$dv <- rec(liss_subj$id, rec = "1=3 [Political Interest]; 2=1 [Vote]; 3=2 [Intention to vote];  
               4=4 [Internal Efficacy]; 5=5 [Pol. part. index]")
liss_subj$term <- as_character(liss_subj$dv)

polat_subj$dv <- rec(polat_subj$id, rec = "3=2 [Political Interest]; 1=1 [Intention to Vote]; 2=7 [Pol. media use]; 4=3 [Internal Efficacy]; 
                5=6 [Duty to vote]; 6=4 [Has party ID]; 7=5 [Party ID strength]; 8=8 [Pol. part. index];  9=9 [Pol. part. onl.]")
polat_subj$term <- as_character(polat_subj$dv)

shp_subj$dv <- rec(shp_subj$id, rec = "1=1 [Vote]; 2=2 [Political Interest]; 3=3 [No. polls part.]")
shp_subj$term <- as_character(shp_subj$dv)

soep_subj$dv <- rec(soep_subj$id, rec = "1=1 [Vote]; 2=3 [Political Interest]; 3=25 [Has party ID]; 
               4=2 [Intention to vote]; 5=26 [Party ID strength]")
soep_subj$term <- as_character(soep_subj$dv)


#us_subj <- bind_rows(gss_subj, nlsy97_subj, psid_subj)
us_subj <- gss_subj
d_subj <- bind_rows(bhps_subj, gles_subj, liss_subj, polat_subj, shp_subj, soep_subj)







#reorder rows for dwplot
bhps_r <- bhps_r[c(9,10,5,6,11,12,7,8,3,4,1,2),]
gles_r <- gles_r[c(21,22,17,18,19,20,15,16,13,14,11,12,9,10,7,8,5,6,3,4,2,1),]
liss_r <- liss_r[c(7,8,5,6,9,10,3,4,1,2),]
polat_r <- polat_r[c(17,18,13,14,7,8,5,6,11,12,15,16,3,4,1,2),]
shp_r <- shp_r[c(5,6,3,4,1,2),]
soep_r <- soep_r[c(9,10,3,4,7,8,5,6,1,2),]
us_r <- us_r[c(9,10,7,8,5,6,3,4,1,2),]

bhps <- bhps[c(9,10,5,6,11,12,7,8,3,4,1,2),]
gles <- gles[c(21,22,17,18,19,20,15,16,13,14,11,12,9,10,7,8,5,6,3,4,2,1),]
liss <- liss[c(7,8,5,6,9,10,3,4,1,2),]
polat <- polat[c(17,18,13,14,7,8,5,6,11,12,15,16,3,4,1,2),]
shp <- shp[c(5,6,3,4,1,2),]
soep <- soep[c(9,10,3,4,7,8,5,6,1,2),]
us <- us[c(9,10,7,8,5,6,3,4,1,2),]

bhps_ln <- bhps_ln[c(9,10,5,6,11,12,7,8,3,4,1,2),]
gles_ln <- gles_ln[c(21,22,17,18,19,20,15,16,13,14,11,12,9,10,7,8,5,6,3,4,2,1),]
liss_ln <- liss_ln[c(7,8,5,6,9,10,3,4,1,2),]
polat_ln <- polat_ln[c(17,18,13,14,7,8,5,6,11,12,15,16,3,4,1,2),]
shp_ln <- shp_ln[c(5,6,3,4,1,2),]
soep_ln <- soep_ln[c(9,10,3,4,7,8,5,6,1,2),]
us_ln <- us_ln[c(9,10,7,8,5,6,3,4,1,2),]

bhps_subj <- bhps_subj[c(9,10,5,6,11,12,7,8,3,4,1,2),]
gles_subj <- gles_subj[c(21,22,17,18,19,20,15,16,13,14,11,12,9,10,7,8,5,6,3,4,2,1),]
liss_subj <- liss_subj[c(7,8,5,6,9,10,3,4,1,2),]
polat_subj <- polat_subj[c(17,18,13,14,7,8,5,6,11,12,15,16,3,4,1,2),]
shp_subj <- shp_subj[c(5,6,3,4,1,2),]
soep_subj <- soep_subj[c(9,10,3,4,7,8,5,6,1,2),]
us_subj <- us_subj[c(3,4,1,2),]


#add labels
bhps_r$model <- rec(bhps_r$model, rec = "1=within; 2=between")
gles_r$model <- rec(gles_r$model, rec = "1=within; 2=between")
liss_r$model <- rec(liss_r$model, rec = "1=within; 2=between")
polat_r$model <- rec(polat_r$model, rec = "1=within; 2=between")
shp_r$model <- rec(shp_r$model, rec = "1=within; 2=between")
soep_r$model <- rec(soep_r$model, rec = "1=within; 2=between")
us_r$model <- rec(us_r$model, rec = "1=within; 2=between")

bhps$model <- rec(bhps$model, rec = "1=within; 2=between")
gles$model <- rec(gles$model, rec = "1=within; 2=between")
liss$model <- rec(liss$model, rec = "1=within; 2=between")
polat$model <- rec(polat$model, rec = "1=within; 2=between")
shp$model <- rec(shp$model, rec = "1=within; 2=between")
soep$model <- rec(soep$model, rec = "1=within; 2=between")
us$model <- rec(us$model, rec = "1=within; 2=between")

bhps_ln$model <- rec(bhps_ln$model, rec = "1=within; 2=between")
gles_ln$model <- rec(gles_ln$model, rec = "1=within; 2=between")
liss_ln$model <- rec(liss_ln$model, rec = "1=within; 2=between")
polat_ln$model <- rec(polat_ln$model, rec = "1=within; 2=between")
shp_ln$model <- rec(shp_ln$model, rec = "1=within; 2=between")
soep_ln$model <- rec(soep_ln$model, rec = "1=within; 2=between")
us_ln$model <- rec(us_ln$model, rec = "1=within; 2=between")

bhps_subj$model <- rec(bhps_subj$model, rec = "1=within; 2=between")
gles_subj$model <- rec(gles_subj$model, rec = "1=within; 2=between")
liss_subj$model <- rec(liss_subj$model, rec = "1=within; 2=between")
polat_subj$model <- rec(polat_subj$model, rec = "1=within; 2=between")
shp_subj$model <- rec(shp_subj$model, rec = "1=within; 2=between")
soep_subj$model <- rec(soep_subj$model, rec = "1=within; 2=between")
us_subj$model <- rec(us_subj$model, rec = "1=within; 2=between")


#multiply effects of LPM by 9 to make estimates comparable
bhps$estimate <- ifelse(bhps$dv == 1, bhps$estimate * 9, bhps$estimate)
bhps$std.error <- ifelse(bhps$dv == 1, bhps$std.error * 9, bhps$std.error)
bhps_subj$estimate <- ifelse(bhps_subj$dv == 1, bhps_subj$estimate * 9, bhps_subj$estimate)
bhps_subj$std.error <- ifelse(bhps_subj$dv == 1, bhps_subj$std.error * 9, bhps_subj$std.error)

gles$estimate <- ifelse(gles$dv == 1 | gles$dv==4, gles$estimate * 9, gles$estimate)
gles$std.error <- ifelse(gles$dv == 1 | gles$dv==4, gles$std.error * 9, gles$std.error)
gles_subj$estimate <- ifelse(gles_subj$dv == 1 | gles_subj$dv==4, gles_subj$estimate * 9, gles_subj$estimate)
gles_subj$std.error <- ifelse(gles_subj$dv == 1 | gles_subj$dv==4, gles_subj$std.error * 9, gles_subj$std.error)

liss$estimate <- ifelse(liss$dv == 1 | liss$dv == 3, liss$estimate * 9, liss$estimate)
liss$std.error <- ifelse(liss$dv == 1 | liss$dv == 3, liss$std.error * 9, liss$std.error)
liss_subj$estimate <- ifelse(liss_subj$dv == 1 | liss_subj$dv == 3, liss_subj$estimate * 9, liss_subj$estimate)
liss_subj$std.error <- ifelse(liss_subj$dv == 1 | liss_subj$dv == 3, liss_subj$std.error * 9, liss_subj$std.error)

polat$estimate <- ifelse(polat$dv == 4, polat$estimate * 9, polat$estimate)
polat$std.error <- ifelse(polat$dv == 4, polat$std.error * 9, polat$std.error)
polat_subj$estimate <- ifelse(polat_subj$dv == 4, polat_subj$estimate * 9, polat_subj$estimate)
polat_subj$std.error <- ifelse(polat_subj$dv == 4, polat_subj$std.error * 9, polat_subj$std.error)

shp$estimate <- ifelse(shp$dv == 1, shp$estimate * 9, shp$estimate)
shp$std.error <- ifelse(shp$dv == 1, shp$std.error * 9, shp$std.error)
shp_subj$estimate <- ifelse(shp_subj$dv == 1, shp_subj$estimate * 9, shp_subj$estimate)
shp_subj$std.error <- ifelse(shp_subj$dv == 1, shp_subj$std.error * 9, shp_subj$std.error)

soep$estimate <- ifelse(soep$dv == 1 | soep$dv == 25, soep$estimate * 9, soep$estimate)
soep$std.error <- ifelse(soep$dv == 1 | soep$dv == 25, soep$std.error * 9, soep$std.error)
soep_subj$estimate <- ifelse(soep_subj$dv == 1 | soep_subj$dv == 25, soep_subj$estimate * 9, soep_subj$estimate)
soep_subj$std.error <- ifelse(soep_subj$dv == 1 | soep_subj$dv == 25, soep_subj$std.error * 9, soep_subj$std.error)

us$estimate <- ifelse(us$dv == 1 | us$dv >= 22, us$estimate * 9, us$estimate)
us$std.error <- ifelse(us$dv == 1 | us$dv >= 22, us$std.error * 9, us$std.error)
#us_subj$estimate <- ifelse(us_subj$dv == 1 | us_subj$dv >= 22, us_subj$estimate * 10, us_subj$estimate)
#us_subj$std.error <- ifelse(us_subj$dv == 1 | us_subj$dv >= 22, us_subj$std.error * 10, us_subj$std.error)

#divide coefficients of loginc models so that it represents 10% changes
bhps_ln$estimate <- ifelse(bhps_ln$dv == 1, bhps_ln$estimate / 10, bhps_ln$estimate / 10)
bhps_ln$std.error <- ifelse(bhps_ln$dv == 1, bhps_ln$std.error / 10, bhps_ln$std.error / 10)

gles_ln$estimate <- ifelse(gles_ln$dv == 1 | gles_ln$dv==4, gles_ln$estimate / 10, gles_ln$estimate / 10)
gles_ln$std.error <- ifelse(gles_ln$dv == 1 | gles_ln$dv==4, gles_ln$std.error / 10, gles_ln$std.error / 10)

liss_ln$estimate <- ifelse(liss_ln$dv == 1 | liss_ln$dv == 3, liss_ln$estimate / 10, liss_ln$estimate / 10)
liss_ln$std.error <- ifelse(liss_ln$dv == 1 | liss_ln$dv == 3, liss_ln$std.error / 10, liss_ln$std.error / 10)

polat_ln$estimate <- ifelse(polat_ln$dv == 4, polat_ln$estimate / 10, polat_ln$estimate / 10)
polat_ln$std.error <- ifelse(polat_ln$dv == 4, polat_ln$std.error / 10, polat_ln$std.error / 10)

shp_ln$estimate <- ifelse(shp_ln$dv == 1, shp_ln$estimate / 10, shp_ln$estimate / 10)
shp_ln$std.error <- ifelse(shp_ln$dv == 1, shp_ln$std.error / 10, shp_ln$std.error / 10)

soep_ln$estimate <- ifelse(soep_ln$dv == 1 | soep_ln$dv == 25, soep_ln$estimate / 10, soep_ln$estimate / 10)
soep_ln$std.error <- ifelse(soep_ln$dv == 1 | soep_ln$dv == 25, soep_ln$std.error / 10, soep_ln$std.error / 10)

us_ln$estimate <- ifelse(us_ln$dv == 1 | us_ln$dv >= 22, us_ln$estimate / 10, us_ln$estimate / 10)
us_ln$std.error <- ifelse(us_ln$dv == 1 | us_ln$dv >= 22, us_ln$std.error / 10, us_ln$std.error / 10)



#### create plots ####

#objective: raw hhinc
p_bhps_r <- dwplot(bhps_r, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.227,0.602), breaks = c(-0.2,0,0.2,0.4,0.6), labels = c("-0.2", "0", "0.2", "0.4", "0.6")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("UK (BHPS & UKHLS)")  
p_bhps_r


p_gles_r <- dwplot(gles_r, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.227,0.602), breaks = c(-0.2,0,0.2,0.4,0.6), labels = c("-0.2", "0", "0.2", "0.4", "0.6")) +  
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Germany (GLES)")

p_liss_r <- dwplot(liss_r, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.227,0.602), breaks = c(-0.2,0,0.2,0.4,0.6), labels = c("-0.2", "0", "0.2", "0.4", "0.6")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Netherlands (LISS)")

p_polat_r <- dwplot(polat_r, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                  whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.227,0.602), breaks = c(-0.2,0,0.2,0.4,0.6), labels = c("-0.2", "0", "0.2", "0.4", "0.6")) +  
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Spain (POLAT)")

p_psid_r <- dwplot(psid_r, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.227,0.602), breaks = c(-0.2,0,0.2,0.4,0.6), labels = c("-0.2", "0", "0.2", "0.4", "0.6")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("PSID")

p_shp_r <- dwplot(shp_r, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.227,0.602), breaks = c(-0.2,0,0.2,0.4,0.6), labels = c("-0.2", "0", "0.2", "0.4", "0.6")) +  
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  guides(col = guide_legend(nrow = 2)) +
  ggtitle("Switzerland (SHP)")

p_soep_r <- dwplot(soep_r, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.227,0.602), breaks = c(-0.2,0,0.2,0.4,0.6), labels = c("-0.2", "0", "0.2", "0.4", "0.6")) +  
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Germany (SOEP)")

p_us_r <- dwplot(us_r, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
               whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.227,0.602), breaks = c(-0.2,0,0.2,0.4,0.6), labels = c("-0.2", "0", "0.2", "0.4", "0.6")) +  
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("United States")

p_legend_r <- dwplot(shp_r, dot_args = list(aes(color = model)), whisker_args = list(aes(color = model)))+
  scale_x_continuous(limits = c(0,0))+
  theme_void()+
  theme(legend.position = c(0.5,0.5),
        legend.key.size = unit(1, "cm"),
        legend.text = element_text(size =  12),
        legend.title = element_blank())
p_legend_r[["layers"]][[1]][["show.legend"]] <- TRUE
p_legend_r


plot_hhinc_r <- ggarrange(p_shp_r, p_liss_r, p_gles_r, p_soep_r, p_polat_r, p_bhps_r, p_us_r, p_legend_r,  nrow = 4, ncol=2, align='v')
ggsave("D:/ScieBo/NWG S�P/Panels/Out/Figures/hybrid_rawinc_col.pdf", plot = plot_hhinc_r, height = 3 , width = 2.5, scale=4.0)
ggsave("D:/ScieBo/NWG S�P/Panels/Out/Figures/hybrid_rawinc_col.emf", plot = plot_hhinc_r, height = 3 , width = 2.5, scale=4.0)




#objective: hhinc_dec_pp_sqrt
p_bhps <- dwplot(bhps, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.11,0.29), breaks = c(-0.1,0,0.1,0.2), labels = c("-0.1", "0", "0.1", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("UK (BHPS & UKHLS)")  
p_bhps


p_gles <- dwplot(gles, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.11,0.29), breaks = c(-0.1,0,0.1,0.2), labels = c("-0.1", "0", "0.1", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Germany (GLES)")
                 
p_liss <- dwplot(liss, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.11,0.29), breaks = c(-0.1,0,0.1,0.2), labels = c("-0.1", "0", "0.1", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Netherlands (LISS)")

p_polat <- dwplot(polat, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                  whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.11,0.29), breaks = c(-0.1,0,0.1,0.2), labels = c("-0.1", "0", "0.1", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Spain (POLAT)")
                 
p_psid <- dwplot(psid, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.11,0.29), breaks = c(-0.1,0,0.1,0.2), labels = c("-0.1", "0", "0.1", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("PSID")
                 
p_shp <- dwplot(shp, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.11,0.29), breaks = c(-0.1,0,0.1,0.2), labels = c("-0.1", "0", "0.1", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  guides(col = guide_legend(nrow = 2)) +
  ggtitle("Switzerland (SHP)")
                
p_soep <- dwplot(soep, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.11,0.29), breaks = c(-0.1,0,0.1,0.2), labels = c("-0.1", "0", "0.1", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Germany (SOEP)")

p_us <- dwplot(us, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
               whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.11,0.29), breaks = c(-0.1,0,0.1,0.2), labels = c("-0.1", "0", "0.1", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("United States")

p_legend <- dwplot(shp, dot_args = list(aes(color = model)), whisker_args = list(aes(color = model)))+
  scale_x_continuous(limits = c(0,0))+
  theme_void()+
  theme(legend.position = c(0.5,0.5),
        legend.key.size = unit(1, "cm"),
        legend.text = element_text(size =  12),
        legend.title = element_blank())
p_legend[["layers"]][[1]][["show.legend"]] <- TRUE
p_legend


plot_hhinc <- ggarrange(p_shp, p_liss, p_gles, p_soep, p_polat, p_bhps, p_us, p_legend,  nrow = 4, ncol=2, align='v')
ggsave("D:/ScieBo/NWG S�P/Panels/Out/Figures/hybrid_eqincdec_col.pdf", plot = plot_hhinc, height = 3 , width = 2.5, scale=4.0)
ggsave("D:/ScieBo/NWG S�P/Panels/Out/Figures/hybrid_eqincdec_col.emf", plot = plot_hhinc, height = 3 , width = 2.5, scale=4.0)



#objective: log income
p_bhps_ln <- dwplot(bhps_ln, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.05,0.116), breaks = c(-0.05,0,0.05,0.1), labels = c("-0.05","0","0.05","0.1")) +    
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("UK (BHPS & UKHLS)")  
p_bhps_ln

p_gles_ln <- dwplot(gles_ln, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                    whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.05,0.116), breaks = c(-0.05,0,0.05,0.1), labels = c("-0.05","0","0.05","0.1")) +    
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Germany (GLES)")

p_liss_ln <- dwplot(liss_ln, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.05,0.116), breaks = c(-0.05,0,0.05,0.1), labels = c("-0.05","0","0.05","0.1")) +      
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Netherlands (LISS)")

p_polat_ln <- dwplot(polat_ln, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                    whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.05,0.116), breaks = c(-0.05,0,0.05,0.1), labels = c("-0.05","0","0.05","0.1")) +       
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Spain (POLAT)")

p_psid_ln <- dwplot(psid_ln, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.05,0.116), breaks = c(-0.05,0,0.05,0.1), labels = c("-0.05","0","0.05","0.1")) +       
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("PSID")

p_shp_ln <- dwplot(shp_ln, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.05,0.116), breaks = c(-0.05,0,0.05,0.1), labels = c("-0.05","0","0.05","0.1")) +       
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  guides(col = guide_legend(nrow = 2)) +
  ggtitle("Switzerland (SHP)")

p_soep_ln <- dwplot(soep_ln, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                 whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.05,0.116), breaks = c(-0.05,0,0.05,0.1), labels = c("-0.05","0","0.05","0.1")) +      
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Germany (SOEP)")

p_us_ln <- dwplot(us_ln, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
               whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.05,0.116), breaks = c(-0.05,0,0.05,0.1), labels = c("-0.05","0","0.05","0.1")) +        
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("United States")

p_legend_ln <- dwplot(shp_ln, dot_args = list(aes(color = model)), whisker_args = list(aes(color = model)))+
  scale_x_continuous(limits = c(0,0))+
  theme_void()+
  theme(legend.position = c(0.5,0.5),
        legend.key.size = unit(1, "cm"),
        legend.text = element_text(size =  12),
        legend.title = element_blank())
p_legend_ln[["layers"]][[1]][["show.legend"]] <- TRUE
p_legend_ln


plot_hhinc_ln <- ggarrange(p_shp_ln, p_liss_ln, p_gles_ln, p_soep_ln, p_polat_ln, p_bhps_ln, p_us_ln, p_legend_ln,  nrow = 4, ncol=2, align='v')
ggsave("./Out/Figures/hybrid_loginc_col.pdf", plot = plot_hhinc_ln, height = 3 , width = 2.5, scale=4.0)
ggsave("./Out/Figures/hybrid_loginc_col.emf", plot = plot_hhinc_ln, height = 3 , width = 2.5, scale=4.0)




#subjective
p_bhps_subj <- dwplot(bhps_subj, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                      whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.1,0.373), breaks = c(-0.1,0,0.1,0.2,0.3), labels = c("-0.1", "0", "0.1", "0.2", "0.3")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("UK (BHPS & UKHLS)") 

p_gles_subj <- dwplot(gles_subj, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                      whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.1,0.373), breaks = c(-0.1,0,0.1,0.2,0.3), labels = c("-0.1", "0", "0.1", "0.2", "0.3")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Germany (GLES)")

p_liss_subj <- dwplot(liss_subj, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                      whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.1,0.373), breaks = c(-0.1,0,0.1,0.2,0.3), labels = c("-0.1", "0", "0.1", "0.2", "0.3")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Netherlands (LISS)")

p_polat_subj <- dwplot(polat_subj, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                       whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.1,0.373), breaks = c(-0.1,0,0.1,0.2,0.3), labels = c("-0.1", "0", "0.1", "0.2", "0.3")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Spain (POLAT)")

p_shp_subj <- dwplot(shp_subj, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                     whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.1,0.373), breaks = c(-0.1,0,0.1,0.2,0.3), labels = c("-0.1", "0", "0.1", "0.2", "0.3")) +  
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Switzerland (SHP)")

p_soep_subj <- dwplot(soep_subj, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                      whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.1,0.373), breaks = c(-0.1,0,0.1,0.2,0.3), labels = c("-0.1", "0", "0.1", "0.2", "0.3")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Germany (SOEP)")

p_us_subj <- dwplot(us_subj, vline = geom_vline(xintercept=0, color = "grey60"), dot_args = list(aes(color = model)),
                    whisker_args = list(aes(color = model))) + 
  scale_x_continuous(limits = c(-0.1,0.373), breaks = c(-0.1,0,0.1,0.2,0.3), labels = c("-0.1", "0", "0.1", "0.2", "0.3")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("United States")


plot_hhinc_subj <- ggarrange(p_shp_subj, p_liss_subj, p_gles_subj, p_soep_subj, p_polat_subj, p_bhps_subj, p_us_subj, p_legend,  nrow = 4, ncol=2, align='v')
ggsave("./Out/Figures/hybrid_satinc_col.pdf", plot = plot_hhinc_subj, height = 3 , width = 2.5, scale=4.0)
ggsave("./Out/Figures/hybrid_satinc_col.emf", plot = plot_hhinc_subj, height = 3 , width = 2.5, scale=4.0)



##### DISTRIBUTED LAG MODELS ####

##### read data ####

#raw hhinc
bhps_lag_r <- read_stata("./BHPS/bhps_coef_lag_hhinc.dta")
liss_lag_r <- read_stata("./LISS/liss_coef_lag_hhinc.dta")
shp_lag_r <- read_stata("./shp/shp_coef_lag_hhinc.dta")
soep_lag_r <- read_stata("./SOEP/soep_coef_lag_hhinc.dta")

#hhinc_dec_pp_sqrt
bhps_lag <- read_stata("./BHPS/bhps_coef_lag_hhinc_dec_pp_sqrt.dta")
liss_lag <- read_stata("./LISS/liss_coef_lag_hhinc_dec_pp_sqrt.dta")
shp_lag <- read_stata("./shp/shp_coef_lag_hhinc_dec_pp_sqrt.dta")
soep_lag <- read_stata("./SOEP/soep_coef_lag_hhinc_dec_pp_sqrt.dta")

#logged income
bhps_lag_ln <- read_stata("./BHPS/bhps_coef_lag_loginc.dta")
liss_lag_ln <- read_stata("./LISS/liss_coef_lag_loginc.dta")
shp_lag_ln <- read_stata("./shp/shp_coef_lag_loginc.dta")
soep_lag_ln <- read_stata("./SOEP/soep_coef_lag_loginc.dta")

#hhinc_sat
bhps_lag_subj <- read_stata("./BHPS/bhps_coef_lag_hhinc_sat.dta")
liss_lag_subj <- read_stata("./LISS/liss_coef_lag_hhinc_sat.dta")
shp_lag_subj <- read_stata("./shp/shp_coef_lag_hhinc_sat.dta")
soep_lag_subj <- read_stata("./SOEP/soep_coef_lag_hhinc_sat.dta")

#incdiff (2 deciles)
bhps_lag_diff <- read_stata("./BHPS/bhps_coef_lag_incdiff2.dta")
liss_lag_diff <- read_stata("./LISS/liss_coef_lag_incdiff2.dta")
shp_lag_diff <- read_stata("./shp/shp_coef_lag_incdiff2.dta")
soep_lag_diff <- read_stata("./SOEP/soep_coef_lag_incdiff2.dta")


####p repare data ####

# delete empty rows
bhps_lag_r <- bhps_lag_r[-11,]
liss_lag_r <- liss_lag_r[-21,]
shp_lag_r <- shp_lag_r[-16,]
soep_lag_r <- soep_lag_r[-16,]

bhps_lag <- bhps_lag[-11,]
liss_lag <- liss_lag[-21,]
shp_lag <- shp_lag[-16,]
soep_lag <- soep_lag[-16,]

bhps_lag_ln <- bhps_lag_ln[-11,]
liss_lag_ln <- liss_lag_ln[-21,]
shp_lag_ln <- shp_lag_ln[-16,]
soep_lag_ln <- soep_lag_ln[-16,]

bhps_lag_subj <- bhps_lag_subj[-11,]
liss_lag_subj <- liss_lag_subj[-21,]
shp_lag_subj <- shp_lag_subj[-16,]
soep_lag_subj <- soep_lag_subj[-16,]

bhps_lag_diff <- bhps_lag_diff[-11,]
liss_lag_diff <- liss_lag_diff[-21,]
shp_lag_diff <- shp_lag_diff[-16,]
soep_lag_diff <- soep_lag_diff[-16,]


#set new colnames 
colnames(bhps_lag_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(liss_lag_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")
colnames(shp_lag_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(soep_lag_r) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")

colnames(bhps_lag) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(liss_lag) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")
colnames(shp_lag) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(soep_lag) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")

colnames(bhps_lag_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(liss_lag_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")
colnames(shp_lag_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(soep_lag_ln) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")

colnames(bhps_lag_subj) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")
colnames(liss_lag_subj) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")
colnames(shp_lag_subj) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(soep_lag_subj) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")

colnames(bhps_lag_diff) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")
colnames(liss_lag_diff) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")
colnames(shp_lag_diff) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(soep_lag_diff) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2",  "id")



#create lag identifiers
bhps_lag_r$lag <- rep(c(0,1,2,3,-1), times=2, each=1)
liss_lag_r$lag <- rep(c(0,1,2,3,-1), times=4, each=1)
shp_lag_r$lag <- rep(c(0,1,2,3,-1), times=3, each=1)
soep_lag_r$lag <- rep(c(0,1,2,3,-1), times=3, each=1)

bhps_lag$lag <- rep(c(0,1,2,3,-1), times=2, each=1)
liss_lag$lag <- rep(c(0,1,2,3,-1), times=4, each=1)
shp_lag$lag <- rep(c(0,1,2,3,-1), times=3, each=1)
soep_lag$lag <- rep(c(0,1,2,3,-1), times=3, each=1)

bhps_lag_ln$lag <- rep(c(0,1,2,3,-1), times=2, each=1)
liss_lag_ln$lag <- rep(c(0,1,2,3,-1), times=4, each=1)
shp_lag_ln$lag <- rep(c(0,1,2,3,-1), times=3, each=1)
soep_lag_ln$lag <- rep(c(0,1,2,3,-1), times=3, each=1)

bhps_lag_subj$lag <- rep(c(0,1,2,3,-1), times=2, each=1)
liss_lag_subj$lag <- rep(c(0,1,2,3,-1), times=4, each=1)
shp_lag_subj$lag <- rep(c(0,1,2,3,-1), times=3, each=1)
soep_lag_subj$lag <- rep(c(0,1,2,3,-1), times=3, each=1)

bhps_lag_diff$lag <- rep(c(0,1,2,3,-1), times=2, each=1)
liss_lag_diff$lag <- rep(c(0,1,2,3,-1), times=4, each=1)
shp_lag_diff$lag <- rep(c(0,1,2,3,-1), times=3, each=1)
soep_lag_diff$lag <- rep(c(0,1,2,3,-1), times=3, each=1)

#add identifiers for datasets
bhps_lag_r$df <- 1
liss_lag_r$df <- 2
shp_lag_r$df <- 3
soep_lag_r$df <- 4

bhps_lag$df <- 1
liss_lag$df <- 2
shp_lag$df <- 3
soep_lag$df <- 4

bhps_lag_ln$df <- 1
liss_lag_ln$df <- 2
shp_lag_ln$df <- 3
soep_lag_ln$df <- 4

bhps_lag_subj$df <- 1
liss_lag_subj$df <- 2
shp_lag_subj$df <- 3
soep_lag_subj$df <- 4

bhps_lag_diff$df <- 1
liss_lag_diff$df <- 2
shp_lag_diff$df <- 3
soep_lag_diff$df <- 4

##add DVs 

#raw inc
bhps_lag_r$dv <- rec(bhps_lag_r$id, rec = "1=1 [Vote]; 2=2 [Political Interest]")
bhps_lag_r$term <- to_factor(bhps_lag_r$dv)
bhps_lag_r$term2 <- as_character(bhps_lag_r$dv)

liss_lag_r$dv <- rec(liss_lag_r$id, rec = "1=2 [Political Interest]; 2=3 [Internal Efficacy]; 3=4 [Pol. Participation];  4=1 [Vote]")
liss_lag_r$term <- to_factor(liss_lag_r$dv)
liss_lag_r$term2 <- as_character(liss_lag_r$dv)

shp_lag_r$dv <- rec(shp_lag_r$id, rec = "1=1 [Vote]; 2=2 [Political Interest]; 3=3 [No. polls part.]")
shp_lag_r$term <- to_factor(shp_lag_r$dv)
shp_lag_r$term2 <- as_character(shp_lag_r$dv)

soep_lag_r$dv <- rec(soep_lag_r$id, rec = "1=1 [Political Interest]; 2=2 [Has party ID]; 3=3 [Party ID strength]")
soep_lag_r$term <- to_factor(soep_lag_r$dv)
soep_lag_r$term2 <- as_character(soep_lag_r$dv)


#inc dec
bhps_lag$dv <- rec(bhps_lag$id, rec = "1=1 [Vote]; 2=2 [Political Interest]")
bhps_lag$term <- to_factor(bhps_lag$dv)
bhps_lag$term2 <- as_character(bhps_lag$dv)

liss_lag$dv <- rec(liss_lag$id, rec = "1=2 [Political Interest]; 2=3 [Internal Efficacy]; 3=4 [Pol. Participation];  4=1 [Vote]")
liss_lag$term <- to_factor(liss_lag$dv)
liss_lag$term2 <- as_character(liss_lag$dv)

shp_lag$dv <- rec(shp_lag$id, rec = "1=1 [Vote]; 2=2 [Political Interest]; 3=3 [No. polls part.]")
shp_lag$term <- to_factor(shp_lag$dv)
shp_lag$term2 <- as_character(shp_lag$dv)

soep_lag$dv <- rec(soep_lag$id, rec = "1=1 [Political Interest]; 2=2 [Has party ID]; 3=3 [Party ID strength]")
soep_lag$term <- to_factor(soep_lag$dv)
soep_lag$term2 <- as_character(soep_lag$dv)

#logged income
bhps_lag_ln$dv <- rec(bhps_lag_ln$id, rec = "1=1 [Vote]; 2=2 [Political Interest]")
bhps_lag_ln$term <- to_factor(bhps_lag_ln$dv)
bhps_lag_ln$term2 <- as_character(bhps_lag_ln$dv)

liss_lag_ln$dv <- rec(liss_lag_ln$id, rec = "1=2 [Political Interest]; 2=3 [Internal Efficacy]; 3=4 [Pol. Participation];  4=1 [Vote]")
liss_lag_ln$term <- to_factor(liss_lag_ln$dv)
liss_lag_ln$term2 <- as_character(liss_lag_ln$dv)

shp_lag_ln$dv <- rec(shp_lag_ln$id, rec = "1=1 [Vote]; 2=2 [Political Interest]; 3=3 [No. polls part.]")
shp_lag_ln$term <- to_factor(shp_lag_ln$dv)
shp_lag_ln$term2 <- as_character(shp_lag_ln$dv)

soep_lag_ln$dv <- rec(soep_lag_ln$id, rec = "1=1 [Political Interest]; 2=2 [Has party ID]; 3=3 [Party ID strength]")
soep_lag_ln$term <- to_factor(soep_lag_ln$dv)
soep_lag_ln$term2 <- as_character(soep_lag_ln$dv)

#subjective inc
bhps_lag_subj$dv <- rec(bhps_lag_subj$id, rec = "1=1 [Vote]; 2=2 [Political Interest]")
bhps_lag_subj$term <- to_factor(bhps_lag_subj$dv)
bhps_lag_subj$term2 <- as_character(bhps_lag_subj$dv)

liss_lag_subj$dv <- rec(liss_lag_subj$id, rec = "1=2 [Political Interest]; 2=3 [Internal Efficacy]; 3=4 [Pol. Participation];  4=1 [Vote]")
liss_lag_subj$term <- to_factor(liss_lag_subj$dv)
liss_lag_subj$term2 <- as_character(liss_lag_subj$dv)

shp_lag_subj$dv <- rec(shp_lag_subj$id, rec = "1=1 [Vote]; 2=2 [Political Interest]; 3=3 [No. polls part.]")
shp_lag_subj$term <- to_factor(shp_lag_subj$dv)
shp_lag_subj$term2 <- as_character(shp_lag_subj$dv)

soep_lag_subj$dv <- rec(soep_lag_subj$id, rec = "1=1 [Political Interest]; 2=2 [Has party ID]; 3=3 [Party ID strength]")
soep_lag_subj$term <- to_factor(soep_lag_subj$dv)
soep_lag_subj$term2 <- as_character(soep_lag_subj$dv)


#shock inc
bhps_lag_diff$dv <- rec(bhps_lag_diff$id, rec = "1=1 [Vote]; 2=2 [Political Interest]")
bhps_lag_diff$term <- to_factor(bhps_lag_diff$dv)
bhps_lag_diff$term2 <- as_character(bhps_lag_diff$dv)

liss_lag_diff$dv <- rec(liss_lag_diff$id, rec = "1=2 [Political Interest]; 2=3 [Internal Efficacy]; 3=4 [Pol. Participation];  4=1 [Vote]")
liss_lag_diff$term <- to_factor(liss_lag_diff$dv)
liss_lag_diff$term2 <- as_character(liss_lag_diff$dv)

shp_lag_diff$dv <- rec(shp_lag_diff$id, rec = "1=1 [Vote]; 2=2 [Political Interest]; 3=3 [No. polls part.]")
shp_lag_diff$term <- to_factor(shp_lag_diff$dv)
shp_lag_diff$term2 <- as_character(shp_lag_diff$dv)

soep_lag_diff$dv <- rec(soep_lag_diff$id, rec = "1=1 [Political Interest]; 2=2 [Has party ID]; 3=3 [Party ID strength]")
soep_lag_diff$term <- to_factor(soep_lag_diff$dv)
soep_lag_diff$term2 <- as_character(soep_lag_diff$dv)



#rescale log income effects so that it refers to a 10% change in income
bhps_lag_ln$estimate <- bhps_lag_ln$estimate  / 10
bhps_lag_ln$std.error <- bhps_lag_ln$std.error / 10
bhps_lag_ln$ci_lower <- bhps_lag_ln$ci_lower / 10
bhps_lag_ln$ci_upper <- bhps_lag_ln$ci_upper / 10

liss_lag_ln$estimate <- liss_lag_ln$estimate / 10
liss_lag_ln$std.error <- liss_lag_ln$std.error / 10
liss_lag_ln$ci_lower <- liss_lag_ln$ci_lower / 10
liss_lag_ln$ci_upper <- liss_lag_ln$ci_upper / 10

shp_lag_ln$estimate <- shp_lag_ln$estimate / 10
shp_lag_ln$std.error <- shp_lag_ln$std.error / 10
shp_lag_ln$ci_lower <- shp_lag_ln$ci_lower / 10
shp_lag_ln$ci_upper <- shp_lag_ln$ci_upper / 10

soep_lag_ln$estimate <- soep_lag_ln$estimate / 10
soep_lag_ln$std.error <- soep_lag_ln$std.error / 10
soep_lag_ln$ci_lower <- soep_lag_ln$ci_lower / 10
soep_lag_ln$ci_upper <- soep_lag_ln$ci_upper / 10









#### LAG PLOTS ####

### raw income

#bhps 
p_bhps_lag_r <- ggplot(bhps_lag_r, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) +
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.17, 0.25), breaks=c(-0.1,0,0.1,0.2), labels=c("-0.1","0","0.1","0.2")) +
  scale_color_discrete(labels = c("Vote", "Political Interest")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) + 
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("UK (BHPS & UKHLS)") 
p_bhps_lag_r


#liss
p_liss_lag_r <- ggplot(liss_lag_r, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) + 
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.17, 0.25), breaks=c(-0.1,0,0.1,0.2), labels=c("-0.1","0","0.1","0.2")) +
  scale_color_discrete(labels = c("Vote", "Pol. Interest", "Int. Efficacy", "Pol. Participation")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  guides(col = guide_legend(nrow = 1)) +
  ggtitle("Netherlands (LISS)") 
p_liss_lag_r



#shp
p_shp_lag_r <- ggplot(shp_lag_r, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) + 
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.17, 0.25), breaks=c(-0.1,0,0.1,0.2), labels=c("-0.1","0","0.1","0.2")) +
  scale_color_discrete(labels = c("Vote", "Political Interest", "No. of Polls Participated")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Switzerland (SHP)") 
p_shp_lag_r


#soep
p_soep_lag_r <- ggplot(soep_lag_r, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) + 
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.17, 0.25), breaks=c(-0.1,0,0.1,0.2), labels=c("-0.1","0","0.1","0.2")) +
  scale_color_discrete(labels = c("Political Interest", "Has Party ID", "Party ID Strength")) +    
  xlab("Time") + 
  ylab("")+ 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Germany (SOEP)") 
p_soep_lag_r




ggarrange(p_bhps_lag_r, p_soep_lag_r, p_liss_lag_r, p_shp_lag_r, align='h', common.legend = F, nrow = 2, ncol=2)

ggsave("./Out/Figures/fe_lags_rawinc_col.pdf", height = 2.5 , width = 3.5, scale=3.2)
ggsave("./Out/Figures/fe_lags_rawinc_col.emf", height = 2.5 , width = 3.5, scale=3.2)




### income deciles

#bhps 
p_bhps_lag <- ggplot(bhps_lag, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) +
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.061, 0.08), breaks=c(-0.04,0,0.04,0.08), labels=c("-0.04","0","0.04","0.08")) +
  scale_color_discrete(labels = c("Vote", "Political Interest")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) + 
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("UK (BHPS & UKHLS)") 
p_bhps_lag


#liss
p_liss_lag <- ggplot(liss_lag, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) + 
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.061, 0.08), breaks=c(-0.04,0,0.04,0.08), labels=c("-0.04","0","0.04","0.08")) +
  scale_color_discrete(labels = c("Vote", "Pol. Interest", "Int. Efficacy", "Pol. Participation")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  guides(col = guide_legend(nrow = 1)) +
  ggtitle("Netherlands (LISS)") 
p_liss_lag


#shp
p_shp_lag <- ggplot(shp_lag, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) + 
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.061, 0.08), breaks=c(-0.04,0,0.04,0.08), labels=c("-0.04","0","0.04","0.08")) +
  scale_color_discrete(labels = c("Vote", "Political Interest", "No. of Polls Participated")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Switzerland (SHP)") 
p_shp_lag


#soep
p_soep_lag <- ggplot(soep_lag, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) + 
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.061, 0.08), breaks=c(-0.04,0,0.04,0.08), labels=c("-0.04","0","0.04","0.08")) +
  scale_color_discrete(labels = c("Political Interest", "Has Party ID", "Party ID Strength")) +    
  xlab("Time") + 
  ylab("")+ 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Germany (SOEP)") 
p_soep_lag




ggarrange(p_bhps_lag, p_soep_lag, p_liss_lag, p_shp_lag, align='h', common.legend = F, nrow = 2, ncol=2)

ggsave("./Out/Figures/fe_lags_eqincdec_col.pdf", height = 2.5 , width = 3.5, scale=3.2)
ggsave("./Out/Figures/fe_lags_eqincdec_col.emf", height = 2.5 , width = 3.5, scale=3.2)



### logged income

#bhps 
p_bhps_lag_ln <- ggplot(bhps_lag_ln, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) +
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.032, 0.05), breaks=c(-0.02,0,0.02, 0.04), labels=c("-0.02","0","0.02", "0.04")) +
  scale_color_discrete(labels = c("Vote", "Political Interest")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) + 
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("UK (BHPS & UKHLS)") 
p_bhps_lag_ln


#liss
p_liss_lag_ln <- ggplot(liss_lag_ln, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) + 
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.032, 0.05), breaks=c(-0.02,0,0.02, 0.04), labels=c("-0.02","0","0.02", "0.04")) +
  scale_color_discrete(labels = c("Vote", "Pol. Interest", "Int. Efficacy", "Pol. Participation")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  guides(col = guide_legend(nrow = 1)) +
  ggtitle("Netherlands (LISS)") 
p_liss_lag_ln


#shp
p_shp_lag_ln <- ggplot(shp_lag_ln, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) + 
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.032, 0.05), breaks=c(-0.02,0,0.02, 0.04), labels=c("-0.02","0","0.02", "0.04")) +
  scale_color_discrete(labels = c("Vote", "Political Interest", "No. of Polls Participated")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Switzerland (SHP)") 
p_shp_lag_ln


#soep
p_soep_lag_ln <- ggplot(soep_lag_ln, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) + 
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.032, 0.05), breaks=c(-0.02,0,0.02, 0.04), labels=c("-0.02","0","0.02", "0.04")) +
  scale_color_discrete(labels = c("Political Interest", "Has Party ID", "Party ID Strength")) +    
  xlab("Time") + 
  ylab("")+ 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Germany (SOEP)") 
p_soep_lag_ln


ggarrange(p_bhps_lag_ln, p_soep_lag_ln, p_liss_lag_ln, p_shp_lag_ln, align='h', common.legend = F, nrow = 2, ncol=2)

ggsave("./Out/Figures/fe_lags_loginc_col.pdf", height = 2.5 , width = 3.5, scale=3.2)
ggsave("./Out/Figures/fe_lags_loginc_col.emf", height = 2.5 , width = 3.5, scale=3.2)




### subjective income

#bhps 
p_bhps_lag_subj <- ggplot(bhps_lag_subj, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) +
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) +  
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.043, 0.04), breaks=c(-0.04,-0.02,0,0.02,0.04), labels=c("-0.04","-0.02","0","0.02","0.04")) +
  scale_color_discrete(labels = c("Vote", "Political Interest")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("UK (BHPS & UKHLS)") 
p_bhps_lag_subj


#liss
p_liss_lag_subj <- ggplot(liss_lag_subj, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) +
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) +  
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.043, 0.04), breaks=c(-0.04,-0.02,0,0.02,0.04), labels=c("-0.04","-0.02","0","0.02","0.04")) +
  scale_color_discrete(labels = c("Vote", "Pol. Interest", "Int. Efficacy", "Pol. Participation")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Netherlands (LISS)") 
p_liss_lag_subj


#shp
p_shp_lag_subj <- ggplot(shp_lag_subj, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) +
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.043, 0.04), breaks=c(-0.04,-0.02,0,0.02,0.04), labels=c("-0.04","-0.02","0","0.02","0.04")) +
  scale_color_discrete(labels = c("Vote", "Political Interest", "No. of Polls Participated")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Switzerland (SHP)") 
p_shp_lag_subj


#soep
p_soep_lag_subj <- ggplot(soep_lag_subj, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) +
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) +  
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.043, 0.04), breaks=c(-0.04,-0.02,0,0.02,0.04), labels=c("-0.04","-0.02","0","0.02","0.04")) +
  scale_color_discrete(labels = c("Political Interest", "Has Party ID", "Party ID Strength")) +    
  xlab("Time") + 
  ylab("")+ 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Germany (SOEP)") 
p_soep_lag_subj




ggarrange(p_bhps_lag_subj, p_soep_lag_subj, p_liss_lag_subj, p_shp_lag_subj,  
          align='h', common.legend = F, nrow = 2, ncol=2)

ggsave("./Out/Figures/fe_lags_satinc_col.pdf", height = 2.5 , width = 3.5, scale=3.2)
ggsave("./Out/Figures/fe_lags_satinc_col.emf", height = 2.5 , width = 3.5, scale=3.2)






### income shock: difference of min. two income deciles

#bhps 
p_bhps_lag_diff <- ggplot(bhps_lag_diff, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) +
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.249, 0.3155), breaks=c(-0.2,-0.1,0,0.1,0.2,0.3), labels=c("-0.2","-0.1","0","0.1","0.2","0.3")) +
  scale_color_discrete(labels = c("Vote", "Political Interest")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("UK (BHPS & UKHLS)") 
p_bhps_lag_diff


#liss
p_liss_lag_diff <- ggplot(liss_lag_diff, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) +
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.249, 0.3155), breaks=c(-0.2,-0.1,0,0.1,0.2,0.3), labels=c("-0.2","-0.1","0","0.1","0.2","0.3")) +
  scale_color_discrete(labels = c("Vote", "Pol. Interest", "Int. Efficacy", "Pol. Participation")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Netherlands (LISS)") 
p_liss_lag_diff


#shp
p_shp_lag_diff <- ggplot(shp_lag_diff, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) +
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.249, 0.3155), breaks=c(-0.2,-0.1,0,0.1,0.2,0.3), labels=c("-0.2","-0.1","0","0.1","0.2","0.3")) +
  scale_color_discrete(labels = c("Vote", "Political Interest", "No. of Polls Participated")) +
  xlab("Time") + 
  ylab("") + 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Switzerland (SHP)") 
p_shp_lag_diff


#soep
p_soep_lag_diff <- ggplot(soep_lag_diff, aes(x=lag, y=estimate, group=term, color=term)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid", color = "grey", size=0.5) +
  geom_point(size=2.1, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin=ci_lower, ymax=ci_upper), linetype=1, size=0.5, width=0, position = position_dodge(width = 0.5)) + 
  theme(legend.position="bottom") +
  scale_y_continuous(limit = c(-0.249, 0.3155), breaks=c(-0.2,-0.1,0,0.1,0.2,0.3), labels=c("-0.2","-0.1","0","0.1","0.2","0.3")) +
  scale_color_discrete(labels = c("Political Interest", "Has Party ID", "Party ID Strength")) +    
  xlab("Time") + 
  ylab("")+ 
  theme(legend.title=element_blank(), legend.margin=margin(0,0,0,0), legend.box.margin=margin(0,0,0,0), plot.title = element_text(size=15), legend.text=element_text(size=12)) +
  guides(shape = guide_legend(override.aes = list(size=2))) +
  ggtitle("Germany (SOEP)") 
p_soep_lag_diff




ggarrange(p_bhps_lag_diff, p_soep_lag_diff, p_liss_lag_diff, p_shp_lag_diff,  
          align='h', common.legend = F, nrow = 2, ncol=2)

ggsave("./Out/Figures/fe_lags_incshock_col.pdf", height = 2.5 , width = 3.5, scale=3.2)
ggsave("./Out/Figures/fe_lags_incshock_col.emf", height = 2.5 , width = 3.5, scale=3.2)




